%  add required paths

addpath('../..')

%  set up a mesh for smoothing weather data over Canada

%  load the lat/long data, lon first in degrees and minutes, and lat next

load LatLongsNum.txt

%  convert to decimal values and reverse latitude

Lng = -(LatLongsNum(:,3) + LatLongsNum(:,4)./60);
Lat =   LatLongsNum(:,1) + LatLongsNum(:,2)./60;

%  load the daily temperature .mat file

load daily

%  plot the points with place names

figure(1)
phdl=plot(Lng,Lat,'o');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{16} Longitude')
ylabel('\fontsize{16} Latitude')
axis([-145, -40, 40, 80])
for i=1:35
    text(Lng(i)+1,Lat(i),place(i,:), 'FontSize', 12)
end

% %  ---------------  Create mesh for observation points ---------------
%  Michelle:  ignore this section, which generates the mesh using 
%  Matlab's code.
% 
% %  Delaunay mesh for all 35 points
% 
% Tri = delaunay(Lng, lat);
% 
% triplot(Tri, Lng, Lat);
% 
% %  drop five stations that are too close to other stations
% 
% %  ptdrop Halifax, Baggotville, Montreal, Toronto, Vancouver
% 
% ptdrop = [2, 9, 12, 14, 26];
% nptdrop = length(ptdrop);
% 
% %  construct the reduced set of points and their indices
% 
% nptleft = 35 - nptdrop;
% LatLngptdrop = zeros(nptleft,2);
% LatLngptindx = zeros(nptleft,1);
% m = 0;
% for i=1:35
%     if ~any(ptdrop == i)
%         m = m + 1;
%         LatLngptdrop(m,1) = Lat(i);
%         LatLngptdrop(m,2) = Lng(i);
%         LatLngptindx(m)   = i;
%     end
% end
% 
% %  replot the mesh with numbered points
% 
% Tridrop  = delaunay(LatLngptdrop);
% nTridrop = size(Tridrop,1);
% 
% figure(1)
% triplot(Tridrop, LatLngptdrop(:,1), LatLngptdrop(:,2));
% for i=1:nptleft
%     Lati = Lat(LatLngptindx(i));
%     Lngi = Lng(LatLngptindx(i));
%     text(Lati+0.5,Lngi,num2str(i));
% end
% 
% figure(2)
% triplot(Tridrop, LatLngptdrop(:,1), LatLngptdrop(:,2));
% for i=1:nptleft
%     Lati = Lat(LatLngptindx(i));
%     Lngi = Lng(LatLngptindx(i));
%     placei = place(LatLngptindx(i),:);
%     text(Lati-1.5,Lngi,num2str(i));
%     text(Lati+0.5,Lngi,placei, 'FontSize', 10)
% end
% 
% %  now drop triangles that cover the exterior of Canada
% 
% for i=1:nTridrop
%     disp([i, sort(Tridrop(i,:)')'])
% end
% 
% tridrop = [6,10,13,30,39,40];
% 
% ntridrop = length(tridrop);
% 
% Tridrop2 = zeros(nTridrop - ntridrop,3);
% 
% m = 0;
% for i=1:ntri
%     if ~any(tridrop == i)
%         m = m + 1;
%         Tridrop2(m,:) = Tri(i,:);
%     end
% end
% 
% %  plot the final mesh with place names
% 
% figure(3)
% triplot(Tridrop2, Lat, Lng);
% 
% for i=LatLngptindx
%     text(Lat(i)+0.5,Lng(i),place(i,:), 'FontSize', 12);
% end
% 
% %  set up the spatial smoothing algorithm
% 
% p = LatLngptdrop;
% e = [];
% t = Tridrop;

%  ---------  Input mesh generated by triangulate in R  -------------------

load Canada_mesh_P.txt
load Canada_mesh_T.txt
load Canada_mesh_E.txt

p = Canada_mesh_P;
e = Canada_mesh_E;
t = Canada_mesh_T;

%  set up the FEM basis object and plot it

order = 1;
% order = 2;

basisobj = create_FEM_basis(p, e, t, order);
nbasis = getnbasis(basisobj);

figure(1)
plot(basisobj)

%  ---------  smooth the data  -------------------

%  set up the monthly precipitaton data

indrng = [  1, 31;
           32, 59;
           60, 90;
           91,120;
          121,151;
          152,181;
          182,212;
          213,243;
          244,273;
          274,304;
          305,334;
          335,365];
      
monthtemp = zeros(35,12);
monthprec = zeros(35,12);
for i=1:12
    monthtemp(:,i) = mean(tempav(indrng(i,1):indrng(i,2),:))';
    monthprec(:,i) = mean(precav(indrng(i,1):indrng(i,2),:))';
end
      
monthlab = ['JAN'; 'FEB'; 'MAR'; 'APR'; 'MAY'; 'JUN'; ...
            'JUL'; 'AUG'; 'SEP'; 'OCT'; 'NOV'; 'DEC'];
        

%  set up the weekly precipitaton data

weekind = zeros(52,7);
for i=1:51
    weekind(i,:) = [(i-1)*7+1 : i*7];
end
weekind(52,1:5) = 51*7+1 : 51*7+5;
weekind(52,6:7) = [1,2];
      
weektemp = zeros(35,52);
weekprec = zeros(35,52);

for i=1:52
    weektemp(:,i) = mean(tempav(weekind(i,:),:))';
    weekprec(:,i) = mean(precav(weekind(i,:),:))';
end
      
monthlab = ['JAN'; 'FEB'; 'MAR'; 'APR'; 'MAY'; 'JUN'; ...
            'JUL'; 'AUG'; 'SEP'; 'OCT'; 'NOV'; 'DEC'];

%  use log10 precipitation as data to be smoothed
%  if precipitation used, surface is dominated by Prince Rupert

weekprecdata = log10(weekprec);
% weekprecdata = weekprec;
weektempdata = weektemp;

%  ---------------  smooth and plot monthy data ---------------------------

%  select the index sequence

% monthindex = 1:12;   %  calendar year
% monthindex = [7:12, 1:6];  % mid-summer to early-summer
monthindex = [5:12, 1:4];  % growing season to early spring
monthindex = 10;

%  smooth the logdata and plot the surface for each month

precwrd = 1;
lambda  = 1e-2;
AZ = 15;
EL = 25;

for imonth=monthindex
    
    if precwrd
%         precdatai = [(1:np)', precdata(LatLngptindx,imonth)];
        precdatai = precdata(:,imonth);
    else
        tempdatai = tempdata(:,imonth);
    end
    
    if order == 2
        lambda = 1e0;
        if precwrd
            [precfdi, laplacefdi] = ...
                smooth_FEM_basis(Lng,Lat,precdatai,basisobj,lambda);
        else
            [tempfdi, laplacefdi] = ...
                smooth_FEM_basis(Lng,Lat,tempdatai,basisobj,lambda);
        end
    else
        if precwrd
            precfdi = smooth_FEM_basis(Lng,Lat,precdatai,basisobj,lambda);
        else
            tempfdi = smooth_FEM_basis(Lng,Lat,tempdatai,basisobj,lambda);
        end
    end
    
    %  plot the smoothed logdata
    
    figure(imonth)
    subplot(1,1,1)
    if precwrd
        plot_FEM(precfdi)
    else
        plot_FEM(tempfdi)
    end
    view(AZ, EL)
    if precwrd
        title(['\fontsize{13} log10 precipitation for month ', ...
               num2str(imonth)])
    else
        title(['\fontsize{13} temperature for month ',num2str(imonth)])
    end
    colorbar
    
%     if order == 2
%         figure(imonth+12)
%         subplot(1,1,1)
%         plot(laplacefdi)
%         title( ...
%         ['\fontsize{13} laplacian of log10 precipitation for month ', ...
%                num2str(imonth)])
%     end
%     colorbar
    
    pause
    
end

%  analyze all the data at once 

%  ---------------------  precipitation  --------------------

lambda  = 1e-2;

[precfd, df, gcv, beta, SSE] = ...
            smooth_FEM_basis(Lng,Lat,precdata,basisobj,lambda);

df
gcv
sqrt(SSE/35)

figure(1)
AZ = 0;
EL = 90;

plot_FEM(precfd, AZ, EL)

monthindex = [5:12, 1:4];  % growing season to early spring
for i=monthindex
    plot_FEM(precfd(i), AZ, EL)
    axis([-150,-50,40,90,-1,1])
    view(AZ,EL)
    caxis([-1,1])
    title(['\fontsize{16} Log_{10} Precipitation for Month ',num2str(i)])
    pause
end

%  ---------------------  temperature  --------------------

lambda = 1e-2;

[tempfd, df, gcv, beta, SSE] = ...
            smooth_FEM_basis(Lng,Lat,tempdata,basisobj,lambda);

df
gcv
sqrt(SSE/35)

figure(2)
AZ = 0;
EL = 90;

plot_FEM(tempfd, AZ, EL, [-35,20], monthlab(monthindex,:))

monthindex = [5:12, 1:4];  % growing season to early spring
for i=monthindex
    plot_FEM(tempfd(i), AZ, EL)
    axis([-150,-50,40,90,-45,25])
    view(AZ,EL)
    caxis([-45,25])
    shading interp
    title(['\fontsize{16} Temperature for Month ',monthlab(i,:)])
    pause
end

%  ---------------  smooth and plot weekly data ---------------------------

%  select the index sequence

% weekindex = 1:5;   %  calendar year
weekindex = [21:52, 1:20];  % growing season to early spring
% weekindex = 10;

%  ------------------------  Temperature  --------------------

lambda = 10^(-0.5);

[weektempfd, df, gcv, beta, SSE] = ...
            smooth_FEM_basis(Lng,Lat,weektempdata,basisobj,lambda);

df
sqrt(mean(gcv))

sqrt(SSE/35)

weektempcoefmat = getcoef(weektempfd);

figure(2)
plot(weekindex, weektempcoefmat(:,weekindex)', 'o')

figure(3)
AZ = 0;
EL = 90;

plot_FEM(weektempfd, AZ, EL, [-45,20])

loglamvec = -3:0.5:2;
nloglam   = length(loglamvec);
dfsave    = zeros(nloglam,1);
gcvsave   = zeros(nloglam,1);

for ilam=1:nloglam
    lambdai = 10^loglamvec(ilam);
    [weektempfdi, dfi, gcvi] = ...
            smooth_FEM_basis(Lng,Lat,weektempdata,basisobj,lambdai);
    dfsave(ilam)  = dfi;
    gcvsave(ilam) = sum(gcvi);
end

[loglamvec',dfsave,gcvsave]

for i=weekindex
    plot_FEM(weektempfd(i), AZ, EL)
    axis([-150,-50,40,90,-45,25])
    view(AZ,EL)
    caxis([-45,25])
    shading interp
    title(['\fontsize{16} Temperature for Week ',num2str(i),' '])
    pause(1)
end

%  Precipitation

lambda = 1;

[weekprecfd, df, gcv, beta, SSE] = ...
            smooth_FEM_basis(Lng,Lat,weekprecdata,basisobj,lambda);

df
sqrt(mean(gcv))

sqrt(mean(SSE))

figure(2)
AZ = 0;
EL = 90;

plot_FEM(weekprecfd, AZ, EL, [-2,2])

loglamvec = -2:0.5:5;
nloglam   = length(loglamvec);
dfsave    = zeros(nloglam,1);
gcvsave   = zeros(nloglam,1);

for ilam=1:nloglam
    lambdai = 10^loglamvec(ilam);
    [weekprecfdi, dfi, gcvi] = ...
            smooth_FEM_basis(Lng,Lat,weekprecdata,basisobj,lambdai);
    dfsave(ilam)  = dfi;
    gcvsave(ilam) = sum(gcvi);
end

[loglamvec',dfsave,gcvsave]

for i=weekindex
    plot_FEM(weekprecfd(i), AZ, EL)
    axis([-150,-50,40,90,-2,2])
    view(AZ,EL)
    caxis([-2,2])
    shading interp
    title(['\fontsize{24} Temperature for Week ',num2str(i),' '])
    pause(1)
end

save WeatherMesh

